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Executive  Summary 


A  semi-analytical  approach  to  three-dimensional  (3-D)  shape  optimization  problems  for  a  viscous 
incompressible  fluid  under  the  assumption  of  zero  (low)  Reynolds  number  has  been  developed.  It 
couples  the  theory  of  generalized  analytic  functions  with  the  adjoint  equations-based  method.  A 
solution  to  Stokes  equations  governing  the  behavior  of  the  fluid  has  been  reduced  to  integral  equations 
based  on  the  Cauchy  integral  formula  for  fc-harmonically  analytic  functions.  The  fluid  velocity  and 
boundary  shape  are  the  state  and  design  variables,  respectively,  and  a  shape  optimization  problem  is  to 
find  shape  minimizing  the  energy  dissipation  rate.  In  contrast  to  the  classical  optimal  control  theory, 
the  shape  optimization  problem  has  been  formulated  as  an  optimal  control  problem  with  constraints 
in  the  form  of  integral  equations.  The  optimality  conditions  (state,  adjoint  and  design  equations) 
for  the  optimal  control  problem  have  been  derived.  The  advantage  of  the  suggested  approach  is 
that  the  state  and  adjoint  variables  are  single-variable  functions,  which  being  represented  analytically 
in  the  form  of  series  with  unknown  coefficients,  can  be  accurately  determined  from  the  state  and 
adjoint  integral  equations,  for  example,  by  minimizing  the  total  squared  error.  The  optimal  shape  has 
been  found  iteratively  by  a  gradient-based  method,  in  which  at  each  iteration,  the  state  and  adjoint 
variables  have  been  determined  for  an  updated  shape  and  the  gradient  for  the  cost  functional  with 
respect  to  the  shape  has  been  obtained  by  the  adjoint  equations-based  method.  The  suggested  semi- 
analytical  approach  has  been  illustrated  for  the  drag  minimization  problem  for  motion  of  a  solid  body 
of  revolution  in  the  viscous  incompressible  fluid  and  has  proved  to  be  efficient  and  accurate. 

The  project  involved  two  Ph.D.  students  Anton  Molyboha  and  Bogdan  Grechuk  from  the  Depart¬ 
ment  of  Mathematical  Sciences,  Stevens  Institute  of  Technology. 
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1  Introduction 


The  theory  of  analytic  functions  of  complex  variable  is  one  of  the  primary  methods  for  analytical 
treatment  of  two-dimensional  (2-D)  problems  in  various  areas  of  applied  mathematics,  including  the 
theory  of  electromagnetism,  elasticity  and  hydrodynamics.  Coupled  with  mathematical  programming 
techniques,  it  provides  an  efficient  framework  for  2-D  optimal  design  problems.  This  project  develops 
the  approach  of  generalized  analytic  functions  in  application  to  3-D  shape  optimization  problems, 
in  particular  for  viscous  incompressible  fluid.  We  demonstrated  the  approach  in  finding  the  optimal 
shape  for  a  solid  unit-volume  body  translating  in  the  fluid  with  the  minimal  resistance  force  (drag). 

The  behavior  of  steady  flows  of  a  viscous  incompressible  fluid  under  the  assumption  of  zero  (low) 
Reynolds  number  (so-called  Stokes’  creeping  flows)  is  described  by  the  Stokes  equations 

H  Au  =  grad  p,  div  u  =  0,  (1) 

where  u  is  the  fluid  velocity  field,  p  is  the  pressure  in  the  fluid,  //  is  the  shear  viscosity,  and  Au  = 
grad(divu)  —  curl  curl  u;  see  [7,  8].  Stokes  flows  about  solid  particles  have  been  and  continue  to  be 
a  popular  subject  of  research  in  fluid  mechanics  [7],  in  particular  in  analytical  hydrodynamics.  The 
drag  and  torque  of  solid  bodies  in  Stokes  flows  are  used  for  designing  chemical  technologies  that 
deal  with  particle  sedimentation;  see  [7].  Special  attention  has  been  paid  to  the  form  of  a  solid 
body  that  minimizes  the  energy  dissipation  rate  in  Stokes  flows.  If  the  body  translates  along  some 
axis,  e.g.,  the  z-axis,  with  a  constant  velocity  then  the  principle  of  minimizing  the  energy  dissipation 
rate  is  equivalent  to  minimizing  the  resisting  force  (drag),  experienced  by  the  body.  In  this  regard, 
the  problem  of  finding  the  optimal  shape  for  the  solid  unit-volume  body  in  the  Stokes  flow  is  a 
benchmark  problem,  which  was  first  addressed  by  Pironneau  [12],  who  established  the  optimality 
condition  ||^||  =  const  and  developed  an  iterative  algorithm  for  finding  the  optimal  shape.  Bourot 
[3]  used  Pironneau’s  algorithm  to  find  the  optimal  shape  and  obtained  the  drag  of  0.95425  compared 
to  that  of  the  unit-volume  sphere.  At  each  iteration,  he  represented  solution  to  Stokes  equations  in 
the  form  of  the  series  of  spheroidal  harmonics1  and  solved  the  boundary  conditions  by  minimizing  the 
total  squared  error  with  respect  to  unknown  coefficients.  However,  Pironneau’s  optimality  condition 
is  not  applicable  for  incorporating  other  constraints  on  the  shape  and/or  motion  of  the  body,  e.g., 
for  translation  of  a  solid  unit-volume  body  of  revolution  in  the  direction  orthogonal  to  its  axis  of 
revolution.  Also,  implementation  of  any  shape  optimization  algorithm  requires  very  accurate  solution 
to  Stokes  equations.  To  solve  this  shape  optimization  problem,  Ogawa  and  Kawahara  [11]  used 
the  finite  element  method  (FEM).  However,  their  result  is  visibly  different  from  the  one  obtained 
by  Bourot.  The  adjoint  equations-based  method  has  proved  to  be  efficient  for  PDE  constrained 
optimization  problems,  in  particular  for  control  problems  in  incompressible  viscous  flows  [6].  It  makes 
use  of  adjoint  equations  to  facilitate  finding  the  gradient  for  objective  function  with  respect  to  design 
variables.  For  additional  details  about  this  shape  optimization  problem  and  adjoint  equations-based 
method,  see  [17,  16,  4,  10,  11,  9,  3,  15,  6,  13,  5,  1,  14], 

The  goal  of  this  work  was  to  develop  efficient  algorithms  for  3-D  shape  optimization  applicable 
to  different  models  of  the  viscous  incompressible  fluid  with  a  variety  of  constraints  on  body’s  shape 
(volume,  surface,  cross-section,  etc.)  and  body’s  motion.  We  developed  the  senri-analytical  approach 
coupling  the  framework  of  generalized  analytic  functions  with  adjoint  equations-based  method  and 
demonstrated  this  approach  for  the  problem  of  finding  optimal  shape  for  the  solid  unit-volume  body 
translating  in  the  fluid. 

In  [20],  we  introduced  a  special  class  of  generalized  analytic  functions  that  arise  from  the  funda¬ 
mental  relationship  between  a  scalar  field  4>  and  vectorial  field  A: 

grad  cj)  =  —  curl  A,  div  A  =  0,  (2) 

1This  solution  form  is  used  for  representing  the  velocity  field  for  prolate  spheroid. 
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which  maintains  that  0  and  A  are  related  scalar  and  vectorial  potentials,  respectively.  This  relationship 
is  encountered  in  various  areas  of  applied  mathematics,  in  particular  in  hydrodynamics,  the  theory  of 
elasticity,  electromagnetism,  etc.;  see  [20].  With  divu  =  0,  the  first  equation  in  (1)  can  be  rewritten 
as  grad  p  =  —//curl  (curlu),  whence  it  follows  that  the  vorticity  u  =  curlu  and  pressure  p  are  related 
by  gradp  =  — // curl u>  with  diva/  =  0  and  thus,  /t a/  and  p  are  related  potentials  satisfying  (2). 

In  2-D  case  in  Cartesian  coordinates,  (2)  reduces  to  the  classical  Cauchy-Riemann  system  for 
ordinary  analytic  functions,  and  in  the  3-D  axially  symmetric  case  in  cylindrical  coordinates  (r,  p,  z )2, 
(2)  defines  so-called  r-analytic  functions;  see  [20].  In  the  3-D  asymmetric  case,  (2)  relates  the  kth 
harmonics  of  <f>  and  A,  k  €  ZjJ",  with  respect  to  p,  and  reduces  to  a  series  of  systems  of  two  linear 
first-order  partial  differential  equations 


A  _  k\  i/(k)  =  d_y(k+ l)  9_  jj{k)  =  _(_d_  +  y(k+ 1) 

dr  r )  dz  ’  dz  \dr  r  ) 


(3) 


which  for  each  k  €  Zq  defines  a  class  of  fc-harnronically  analytic  functions  G^k\r,z)  =  U^k\r,z)  + 
i y(k+1)(r,z),  where  i  =  \/^l;  see  [20]. 

It  follows  from  (3)  that  U ^  and  V('k+1')  are  /c-harmonic  and  (£;+l)-harmonic  functions,  respectively, 
i.e. ,  they  satisfy 

A  kU{k)  =  0  and  Ak+1V{k+1)  =  0,  (4) 

where  A&  denotes  the  so-called  /c-harmonic  operator: 

_  d2  1  d  d2  k2  . 

k  dr 2  r  dr  dz2  r2 

In  [18,  19],  we  developed  a  unifying  framework  for  Ai-harmonically  analytic  functions  in  application 
to  3-D  Stokes  flow  problems  (the  papers  have  been  submitted  for  publication).  In  particular,  we 
obtained  representations  for  the  fluid  velocity  field  in  terms  of  /c-harnronically  analytic  functions 
in  axially  symmetric  and  asymmetric  cases  and  also  derived  a  generalized  Cauchy  integral  formula 
for  A:-harmonically  analytic  functions,  which  paved  the  way  for  reducing  3-D  Stokes  flow  problems 
to  integral  equations  for  /c-harmonically  analytic  functions.3  We  formulated  the  benchmark  shape 
optimization  problem  as  an  integral  equation  constrained  optimization  problem  and  solved  it  using 
the  adjoint  equations-based  method.  As  another  illustration,  we  applied  the  developed  approach  to 
shape  optimization  problem  for  the  solid  unit-volume  body  of  revolution  translating  in  the  fluid  in 
the  direction  perpendicular  to  its  axis  of  revolution.  The  novelty  of  the  suggested  approach  is  in  its 
accuracy  and  efficiency. 


2  Problem  Formulation  and  Optimality  Conditions 

Let  a  solid  unit-volume  body  translate  in  a  viscous  incompressible  fluid  under  the  assumption  of  zero 
(low)  Reynolds  number  with  constant  velocity  vz  along  the  z-axis.  In  this  case,  the  fluid  velocity  field 
is  governed  by  the  Stokes  equations  (1)  and  satisfies  the  no-slip  boundary  conditions  on  the  surface  S 
of  the  body  and  conditions  at  infinity: 

u\s  =  vzk,  u|oo  =  0,  Ploo  =  0.  (6) 

2In  this  case,  the  2-axis  is  the  axis  of  symmetry,  and  <f>  and  A  are  independent  of  the  angular  coordinate  ip. 
3Another  possible  approach  to  solve  boundary-value  problems  for  generalized  analytic  functions  is  to  use  formal 
powers  in  the  sense  of  Bers  [2]. 
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Let  be  the  domain  exterior  to  the  body  T>+  (S  =  dT>+).  The  benchmark  shape  optimiza¬ 
tion  problem  is  to  minimize  the  energy  dissipation  rate  in  the  domain  D_  subject  to  the  boundary 
conditions  (6)  (PDE  constrained  optimization): 


min 

av _ 

JJJ  (curlu) 2  dV 

s.t. 

Stokes  equations  (1), 
boundary  conditions  (6), 

(7) 

unit  volume:  JJJ  dV  =  1. 

Since  the  body  translates  with  the  constant  velocity  we  can  represent  the  energy  dissipation 
rate  via  the  resisting  force  (drag)  F  in  the  direction  of  the  2-axis 


(curlu)2  dV  =  —  uzk  •  F  =  —  vz  Fz 


V- 


It  follows  from  the  symmetry  of  Fz,  Stokes  equations  (1)  and  boundary  conditions  (6)  that  the 
2-axis  is  body’s  axis  of  revolution  and  thus,  the  fluid  velocity  field  is  axially  symmetric,  see  Figure  1. 


<P 


|v,k 


Figure  1:  Translation  of  the  solid  unit- volume  body  in  the  fluid  with  the  velocity  vz  along  the  2-axis. 

Let  ( r,(p,z )  be  the  cylindrical  system  of  coordinates  with  the  basis  (e^e^k)  and  let  1+  be  the 
boundary  of  the  body  T>+  in  the  r 2-cross-sectional  plane  for  r  >  0.  The  components  ur  and  uz 
of  the  axially  symmetric  velocity  field  u  are  independent  of  the  angular  coordinate  <p  and  =  0. 
Consequently,  u(r,  ip,  z )  =  ur(r,  z )  er  +  uz(r,  z)  k.  We  also  introduce  the  complex  variable  (  =  r  +  iz. 
In  [18],  we  represented  the  axially  symmetric  velocity  field  u  in  terms  of  two  0-harmonically  analytic 
functions  G*i(C)  and  G^C)4: 


uz  +  iur  =  (z  -  i  r)  Gi(C)  +  G2(Q.  (8) 

Thus,  the  shape  optimization  problem  is  to  find  the  boundary  such  that  minimizes  the  drag5 

4In  our  context,  the  notation  G(Q  does  not  assume  analyticity  of  G  and  simply  denotes  G(r,  z). 

5In  [18],  we  expressed  the  drag  in  terms  of  the  function  G\. 
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subject  to  the  boundary  conditions  formulated  in  terms  of  the  functions  G\  and  G 2- 


min  Re  |  —  J  ?’Gi(C)dc| 


s.t.  Cri(C)  and  G2(C)  are  O-harmonically  analytic  functions  in  X>_, 

(z-  lr)G1(0  +  G2(0=vz,  (€£+, 


(9) 


unit  volume:  27r  //  r  dr  dz  =  1. 

u+ 


There  are  two  issues  in  dealing  with  (9):  (i)  satisfying  the  condition  that  G'i(C)  and  G2(Q  are 
0-harnronically  analytic  functions  in  D_,  and  (ii)  finding  unknown  shape  l+. 

The  boundary  values  of  Gi(£)  and  G2(C)  at  £+  should  satisfy  the  Sokhotsky-Plenrelj’s  formula, 
which  follows  from  the  Cauchy  integral  formula  for  /c-harnronically  analytic  functions  derived  in  [18]: 


(2  -  ^  )  GAO  +  ±1  Gj(r)  iT  _  a Ar) 


Q^(C  t) 

-2'  df  =  0,  (e£+,  l<j<2,  (10) 

r  +  C 


where  and  f2^(C,  t)  are  real-valued  functions  (see  [18])  and  cc(C)  is  the  angle  between  the 

right  and  left  tangent  vectors  to  the  curve  i+  at  the  point  (.  For  all  points  in  which  £+  is  smooth, 
a(()  =  7 r. 

In  our  work  [18],  we  reduced  the  boundary-value  problem  (2  —  |  r)  Gi(C)  +  G2(()  =  Vz  011  f°r 
two  0-harnronically  analytic  functions  Gri(£)  and  G2(£)  in  the  domain  D _  to  a  single  integral  equation: 


C 


(Gi)  =  jT  (Gi(t)  Ki{Ct)  dT  -  Gi(r)  A2(C,t)  dr)  =  2uz,  C  G  ^ 


where  £  =  r  +  i  z,  r  =  rq  +  i  zi  and 


Ai(C,r)  =  ^  C-xCC,  r)nf(C,r),  tf2(C,r)  =  A  C2(C,  r)  0L°J(C,  r) 


(0)/ 


(11) 


(12) 


with 


Gi(C,r)  = 


zi  -  ±n-(z-  Jr) 


T~  C 

Thus,  the  problem  (9)  reduces  to 


C2(C,t)  = 


zi+  \  rx-{z-  \  r) 

C  +  t 


min  Re  J  ?’Gi(C)dc| 
s.t.  C(Gi)  =  2vz,  (  G  £+, 
27 r  U  r  dr  dz  =  1. 


(13) 


Let  the  curve  ■£+  be  parameterized  by  s  E  [0, 1],  i.e.,  £(s)  =  r(s)  +  i^(s),  and  let  G\  be  determined 
as  a  function  of  s,  i.e.,  Gi  =  G*i(s).  Then  the  problem  (13)  is  reformulated  as  an  optimal  control 
problem: 


min  Re  <  —  /  r(s)  Gi(s)  C(s)  ds 

COO  l  Jo 

s.t.  S^(GiC)  =  2vz,  s  G  [0,1], 
7T  /  r2(s)  z(s)  ds  =  1, 

Jo 

r(0)  =  0,  r(l)  =  0, 


(14) 
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where  the  first  constraint  is  the  state  equation  and  the  operator  is  determined  by 

rl 


SdGi  0  = 

We  also  assume  that 


'o 


(15) 


lim  (r(s)  Gi(s))  =  0,  lim  (r(s)  G\(s))  =  0. 

s— >1 


Let  A(s)  be  a  complex-valued  function  and  let  r/  be  a  real-valued  constant.  Then,  the  Lagrangian 
for  (14)  is  determined  by 


L{(,  C;  G\  ;  A,  rj)  =  Re 
where 


-  r(s))Gi(s)C(s)  -  2vz  A (s)  +  rj r2(s)  z(s) )  ds  - 


A(A)  =  f  (mKiimXis))  -  \(t)K2(((t)X(s)))  dt. 


In  this  case,  the  adjoint  (or  co-state)  equation  takes  the  form 

v4c(A)=r(s)|  s€[0,l]. 

With  (18),  the  total  variation  of  the  Lagrangian  (16)  reduces  to 


(16) 


(17) 


(18) 


5L  =  Re 


1  Gi(s)  C(s)  5^C(A)  r(s))  +  A(s)  dS^G iQ  +  2r?r(s)  ]  5r(s)  ds 


'  o 


dr 


dr 


+  f[  Gi(s)  C(s)  a(,4c(  r(g))  +  A(s)  95clGl°  -  2r]r(s)  r(s)  )  Sz(s)  ds 


(19) 


/  o 


dz 


Optimality  conditions  are  summarized  below: 

rl 


State  eq.: 


Adjoint  eq.: 


/o 


*i(CW,C(t))  -  G!(t)C(i)  A'2(C(^),C(t))  eft  =  2u„  (20a) 


5c(GiC) 


~  X{t)K2{C{t),C{s))^  dt  =  r(s),  (20b) 


Ac  (A) 

Design  eq.’s:  Re  { Gi(s)  ( (s)  +  A(s)  +  2rjr(s)  i(s)]-  =  0,  (20c) 


dr 


Re  J  Gds)  C(s)  ^c(Aj  r(5))  +  A(a)  a<SciGl^  -  2r)  r(s)  r(s)  \  =  0,  (20d) 


Boundary  cond.’s: 


dz  dz 

7 r  f  r2(s )  z(s)  ds  =  1, 

Jo 

r(  0)  =  0,  r(l)  =  0. 


(20e) 

(20f) 
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Observe  that  equations  (20a)-(20d)  are  dependent: 


,d(Ac(X)-r(s))  ,  w  ^(GbC)  ,  0  ,  w  , 

r{s)  Gi(s)  C (s) - ^ - b  A (s)  — Ar - b  2 rjr(s)  z{s ) 


dr 


dr 


,  L  ,^MW)-r(s))  ,  w  .dS^GiQ  0  ,  w, 

+  z(s)  |  Gi(s)  C (s) ^ - b  A(s)  — ^ - 2rjr(s)  r{s ) 


<9z 


/~i  t  \  }- 1  ,d(Ac(X)-r(s))  ,  w  ,dSc(GiC)  n 
=  Gi(s)  C(s) - - — j; - +  A(s)  — ^ - =  0, 


ds 


ds 


whence  it  follows  that 


^  ^  5(^(A)  -  r(s))  \ {  \  dS^GA)  o  /  \  •  t  \ 

Gi{s)  C(s) - ^ - b  A(s) — - 2 nr{s)  r{s) 


r{s) 

z(s) 


dz 

,d(Ac(X)~r(s))  ,  w  ,3Sc(GiC)  ,  n  ,  w 
Gi(s)  C(s) - h - b  A(a)  — ^ - b  2 nr(s)  z(s)  . 


dr 


dr 


Solving  the  system  (20a)-(20f)  analytically  is  still  an  open  issue.  However,  the  optimization  prob¬ 
lem  (14)  can  be  efficiently  solved  by  the  adjoint  equations-based  method,  which  is  the  subject  of  the 
next  section. 


3  Adjoint  Equations-based  Method 

The  problem  (14)  can  be  rewritten  in  the  following  form 

min  Re|  — (GiC,r)| 

COO  ^  ) 

s.t.  S^(GiQ  =  2vz,  s  G  [0, 1], 

C(s)  £  bf , 

where  the  inner  product  (•,  •)  is  in  £2([0, 1]),  i.e. ,  is  defined  by 

(f,g)=[  f(t)g(t)dt ,  ||/||  =  VifJ) 

Jo 


and 


Ab  =  \  ( r(s),z(s )) 


7 tJ  r2(s)z(s)ds  =  1,  r(0)  =  0,  r(l)  =  o| 


(21) 


(22) 


Then  the  Lagrangian  for  (21)  is  given  by 

L(C,CGi,  A)  =  Re{  — (GbC,r)  +  (5C(GX C)  -  2u„  A)} 

=  Re  { (Gi  C,  <S<?(A)  —  r)  —  (2uz,  A) }  , 

where  £(s)  G  X  and  is  the  operator  adjoint  to  S ^  and  is  determined  by  <S£( A) 
Appendix  A  for  derivation  of  the  adjoint  operator). 

The  total  variation  of  the  Lagrangian  takes  the  form 

dL(C,CGi,A)  =Re{<GiC,dc(Sc*(A)  -  r))  +  (6(G1  Q,S*C(X)  -  r) 

+  (S^GiQ  —  2vz,  dA}},  £(s)eA, 


(23) 
A  (A)  (see 


which  by  the  adjoint  equations-based  method  reduces  to  the  variation  only  with  respect  to  (: 


SLiCCG,,  A)  =  S(L  =  Reft^GiC),  A>  -  (G^^r)} 


where 


S&V)  = 


s.t.  state  eq.:  S^(G\()  =  2vz, 

adjoint  eq.:  5^(A)  =  r, 

C(s )  £  X i 


SKdCisMW-m  SK2(C(s), 


STr  /A  \  9Kj( C,T)  dKj(C,T)  dKj(C,T)  dKj((,T)  . 

5Kj(C,T)  =  — J- - Sr-\ - - Sri  H - w - Sz  H - ^ - 8z i,  j  =  1,2. 

or  ori  02:  ozi 

The  corresponding  derivatives  of  K\  and  I\ 2  are  presented  in  Appendix  B. 

The  variation  of  the  Lagrangian  is  used  to  find  the  direction  SC  in  which  the  drag  value  has  the 

steepest  descent.  Let 

F(5C)  =  Re{(^c5c(G1C),A>  -  <GiC»}, 

which  is  a  linear  operator  with  respect  to  S(,  and  let 


C(s)  =  7C(s)  =  7yfe  aj  5(a)  €  X,  (24) 

where  Cj(s)  are  basis  functions  and  7  =  7(01, . . . ,  an )  is  the  multiplier  chosen  to  satisfy  the  unit-volume 
constraint.  Then  we  can  write 

S(  =  Y^j= +  =  S( s )> 

where  bj  =  a,j  ^7=1  +  7 4  a-/,  or  in  vectorial  notation,  if  a  =  (oq, . . . ,  an)T,  d  =  (5ai, . . . ,  San)T 

and  b  =  (61, ... ,  bn)T ,  then  b  =  (V7  •  d)  a  +  7d.  In  terms  of  (24),  the  unit- volume  constraint  takes 
the  form  7T73  r2(s)  z(s)  ds  =  1  and  thus,  the  gradient  of  7  with  respect  to  a  is  determined  by 


<97  Jo  (2  fj(s)f(s)z(s)  +  f2(s)zj(s 

dcij  n  1  /q  (  rl  ^0/  \  \  7  \  ^3 


37T1/3  (  L  7^(5)  z(s)  ds 


From  (25),  we  have  (V7  •  a)  =  —7,  whence  (V7  •  b)  =  (V7  •  d)(V7  •  a)  +  7(V7  •  d)  =  0. 

To  find  the  direction  d*  of  the  steepest  descent,  we  formulate  the  following  problem 

min  F(5C)  =  nrin  (f  •  b)  subject  to  b  7Lb  =  l,  (V7  •  b)  =  0,  (26) 

where  the  vector  f  has  components  fj  =  F(Cj ),  H  is  the  symmetric  matrix  with  elements  H{j  =  (Q,  Cj) 
and  Gi(s)  and  A(s)  are  found  from  the  state  and  adjoint  equations,  respectively  (this  will  be  discussed 
in  the  end  of  the  section).  Solving  (26),  we  obtain 


—  1IU2 \  “I 


\7^TH  1V7  J  \\7^tH  1V7 


g-‘V7  -  H-‘f)  , 


and  consequently,  we  can  express  d  via  b  as 


d*  =  -(b*  -  (V7  •  d*)  a) . 
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Since  the  scalar  product  of  (27)  with  V7  reduces  to  identity,  d*  can  be  chosen  so  that  (V7  •  d*)  =  0. 
Then  the  new  vector  anew  is  determined  by 

&new  —  a  "b  h  d  , 

where  the  step  size  h  is  found  by  the  golden  section  technique. 

Finally,  representing  the  functions  Gi(s)  and  A(s)  on  [0, 1]  in  the  form 

Era  ^  ^ 

,  Apikdlk(s)  +  Xp2k92k(s)), 

^  ~  ~  (28) 

A(s)  =  />,  Aqik^ik(s)  +  \q2k  A2fc(s)), 

z — *k=  1 

where  gik(s),  g2k(s ),  Aifc(s)  and  A2 k{s)  are  basis  functions,  satisfying  the  boundary  conditions  r(0)  =  0 
and  r(l)  =  0  and  the  symmetry  conditions  Gi(l  — s)  =  — Gi(s)  and  A(l  — s)  =  —  A(s),  we  find  unknown 
coefficients  p\k_,  P2k->  Qik  and  q2k  by  minimizing  the  total  squared  error 

min  ||5c(GiC)  —  2vz\\2,  min  ||5?(A)  —  r||2,  (29) 

pik,P2k  qik,q2k  ^ 

where  ||  •  ||  is  the  norm  in  £2([0, 1]). 

Let  the  coefficients  pn,. . .,  pim,  P21,  ■  ■  P2m  form  single  vector  p  and  let  a  matrix  J  consist  of 
four  blocks,  each  of  which  having  the  components 

Jjk  =  Re{(5c(C?ij(s)),5c(C?ifc(s)))},  J-i  =  Re{(5c(C5ij(s)),5c(iC?2it(s)))}, 

Jjk  =  Jlji  Jjk  =  M(‘5c(iC?2i(s)),5c(iC?2fc(s)))}- 

Also  let  components 

wk  =  2u2Re{(l,5c(C?ifc(s)))},  wm+k  =  2vz  Re{(l,  5c(iC?2fe(s)))},  1  <  k  <  m, 

form  vector  w,  then  the  first  problem  in  (29)  reduces  to  a  simple  quadratic  optimization  problem 

min  pT  Jp  -  2w  p,  (30) 

p 

whose  solution  is  given  by  p  =  J-1  w.  The  second  minimization  problem  in  (29)  is  solved  similarly. 
We  summarize  the  adjoint  equations-based  method  for  finding  the  optimal  shape. 

Algorithm  1 

1.  Parametrize  (( s )  =  7(01, . . . ,  an )  J2j= 1  aj  Cj(s)  £  A  and  set  initial  a  =  (ai, . . . ,  an), 

2.  Solve  the  state  and  adjoint  equations  for  a  given  shape  (  by  minimizing  the  total  squared  error 
(problem  (29) ), 

3.  Calculate  the  optimal  direction  d*  by  (27), 

4-  Find  the  optimal  step  h  in  the  direction  d*  and  update  a, 

5.  If  1 1  hd*  1 1  <  e  then  Stop,  otherwise  Go  To  Step  2. 

The  next  section  implements  the  algorithm  and  presents  computational  results. 
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4  Computational  Results 

Let  body’s  shape  be  represented  in  the  form 

ETt  ■  

k=iake-msT2k-2(2s-l),  a  €[0,1],  (31) 

where  Tfc(i)  are  Chebyshev’s  polynomials  of  the  first  kind  and  the  multiplier  7  =  7(01, . . . ,  an)  is  found 
from  the  constraint  7T73  z(s )  ds  =  1  (unit  volume).  Obviously,  r(0)  =  0  and  r(l)  =  0. 

Solutions  to  the  state  and  adjoint  equations  are  found  from  (28)  with 

9ik(s)  =  T2k-i(2s  -  1),  ?2fc(s)  =  T2k-2{2s  -  1),  1  <  k  <  m, 

Xik(s)  =  T2fc_i(2s  -  1),  hk(s)  =  T2fc_2(2s  -  1),  1  <  k  <  m. 

We  solved  the  shape  optimization  problem  for  m  =  12  starting  from  the  unit-volume  sphere.  For 
the  first  iteration,  we  used  n  =  2  in  the  basis  {Ck(t)Yk=i  and  then  after  each  3  iterations,  we  increased 
n  by  1.  This  adaptive  basis  procedure  proved  to  be  very  efficient.  We  obtained  the  drag  value  of 
0.95426  after  8  iterations  (this  value  is  normalized  to  the  drag  of  the  unit-volume  sphere).6  Figure 
2  plots  the  objective  function  against  iteration  number.  For  example,  after  13  iterations  (n  =  6),  we 
obtained  the  drag  of  0.954258  and  the  following  shape: 

7  =  0.5988558320119997, 

01  =  1.294644111408456,  a2  =  0.4790126176431742,  a3  =  0.02076294285722642, 

a4  =  -0.008671329483245576,  a5  =  -0.002115187983386224,  a5  =  -0.00024273396227855588. 

Figure  3  shows  the  optimal  shape  for  the  solid  unit- volume  body  translating  in  the  fluid  with  constant 
velocity  and  minimal  drag.  The  shape  is  (almost)  identical  to  the  one  obtained  by  Bourot  [3]  and  is 
a  significant  improvement  compared  to  the  result  of  Ogawa  and  Kawahara  [11]. 

As  expected,  the  performance  of  the  algorithm,  which  is  gradient-based,  depends  on  the  initial 
choice  of  the  shape  and  on  the  method  for  finding  optimal  step  size  for  a  given  direction.  As  an  im¬ 
provement  for  the  algorithm,  we  can  employ  a  conjugate  gradient  method,  which  uses  the  information 
about  several  consecutive  gradients. 


Figure  2:  Value  of  the  objective  function,  normalized  to  the  drag  of  the  unit-volume  sphere,  for  each 
iteration.  The  drag  value  of  0.95426  was  obtained  after  8  iterations  (rn  =  12  and  n  =  5)  by  adaptive 
basis  procedure. 


GBourot’s  value  is  0.95425;  see  [3]. 
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Figure  3:  Optimal  shape  for  the  solid  unit-volume  body  translating  in  the  fluid  with  constant  velocity 
and  the  drag  of  0.95426  (rri  =  12,  n  =  5,  8  iterations). 


5  Optimal  Shape  for  Transversal  Translation 

In  this  section,  we  solve  the  problem  of  finding  optimal  shape  for  a  solid  unit-volume  body  of  revolution 
translating  in  a  viscous  incompressible  fluid  in  the  direction  transversal  to  the  axis  of  revolution.  To 
the  best  of  our  knowledge,  this  problem  has  not  been  addressed.  As  in  the  previous  case,  we  minimize 
the  energy  dissipation  rate  subject  to  the  condition  that  the  fluid  velocity  field  is  governed  by  the 
Stokes  equations  (1).  Here,  we  assume  that  the  z-axis  is  body’s  axis  of  revolution  and  that  the  body 
translates  along  the  y-axis  with  the  constant  velocity  vyy,  see  Figure  4.  In  this  case,  the  boundary 
conditions  for  the  velocity  field  are  given  by 

uls  =  ^j,  u|oo  =  0,  p|oo  =  0.  (33) 


Figure  4:  Transversal  translation  of  the  solid  unit-volume  body  of  revolution  in  the  fluid. 

The  shape  optimization  problem  is  formulated  similarly  to  the  problem  (7)  (PDE  constrained 
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optimization): 


min 

&D- 

s.t. 


(curlu)2  dV 


JJJv _ 

Stokes  equations  (1), 
boundary  conditions  (33), 


unit  volume:  jjj  dV  =  1, 
z-axis  is  body’s  axis  of  revolution, 


where  T>+  and  denote  interior  and  exterior  domains  for  the  body,  respectively. 

As  in  the  previous  case,  since  the  body  translates  with  the  constant  velocity  vy  j ,  we  can  represent 
the  energy  dissipation  rate  via  the  resisting  force  (drag)  F  in  the  direction  of  the  y-axis 


(curlu)2  dV 


-Vy3  F  =  ~vy  Fy 


In  our  work  [19],  we  reduced  the  Stokes  equations  (1)  for  the  boundary  conditions  (33)  to  two 
integral  equations  using  Sokhotsky-Plemelj’s  formula  (10)  for  /c-harmonically  analytic  functions.  Con¬ 
sequently,  the  shape  optimization  problem  reduces  to  minimizing  the  drag  subject  to  two  integral 
equation  constraints: 

min  Re  (  — (Gi0)  C,  r)} 

COO  >-  J 

s.t.  G^\s)  =  H1(s)  +  iH2(s) 

G2(s)  =  \rH2(s)  +iif3(s),  (34) 

TdG{? 0  +  MG2O  =  -4vy,  s  €  [0, 1], 

Re{/3(s)  G2(s)  +  TZ^(G2Q}  =  0,  s£  [0, 1], 

C(s)  €  X, 

where  C  =  r +i  z  and  r  =  r\  +i  z\  are  complex  variables  in  the  rz-cross-sectional  plane  in  the  cylindrical 
coordinates  (r,  <p,  z ),  in  which  the  z-axis  coincides  with  the  z-axis  in  the  Cartesian  coordinates  (x,  y,  z) ; 
Hi(s),  H2(s)  and  H3(s)  are  unknown  real-valued  functions7  parametrized  by  s  on  [0,1];  P(s)  = 
2  —  +C+))  with  a(C)  same  as  in  (10);  and 

T?(Gfc)  =  £  (Gf\t)dt)M1{as),at))  -  dt , 


VdGoX)  = 


dt, 


Af1(C,r)  =  -C11(C,r)^0)(C,r), 

7Tt 


Cn(C) r)  = 


z  —  z\  —  2  i{r  —  ri) 

t-C, 


7T1 


T~( 


M2(C,t)  =  -C'21(C,t)^°)(C,t), 

7TI 

z  —  z\  —  2  i(r  +  ri) 


C'2i(C,t)  = 


r  +  C 


,  1  n^\c,r)  +  n^\c,T) 

M4(C,t)  =  — 


7Tt 


r  +  C 


7 In  [19],  the  functions  H2G)  and  fAjs)  correspond  to  U[°\s),  R/^(s)  and  Ug°\s),  respectively. 
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V(< 

W  0  V 


KdG2C)  =  -  G2WCW 


n$\as),m 
m  -  as) 


+  G2(t)  a t ) 


ng(CM,CW)' 
CW  +  C(S) 


dt, 


m  Jo  \ 

and  the  set  X  is  defined  by  (22). 

We  solved  the  problem  (34)  using  the  adjoint  equations-based  method  discussed  in  Section  3. 

Let  Ai(s)  £  C  and  A2(s)  £  M,  then  the  Lagrangian  for  (34)  takes  the  form 

L( a  H3,Ti,  A2)  =  Re  {-(gS0)  C,r>  +  (Tc(Gf  }C)  +  VC(G2 C)  -  4 vy,Ti)  +  (PG2  +  ftc(G2C),  A2)} 

=  Re{(F!,C(^(Ai)  -r))  +  (H3,  -i?(^(Ai)  +  7^(A2))  -  i/3  A2)  -  4%(Ai,  1) 
+(if2,-iC(Tc*(Ai)-r)+  |C(^(Ai)+^(A2))+  §, 

whence  it  follows  that  the  adjoint  equations  are 


Re 


{C(T?*(A!)— r)}=0, 


Re{-iC(^c(A1)-r)+  5C(^c(Ai)+^c(a2))+  § 

Re  {-i^(Ai)  +  7^(A2))  -  i/3  A2}  =  0, 


=  0, 


or  equivalently, 


Re 


{cc?  (Ar)  -r)}  =  0,  C  (rc*(Ai)  -  r  -  ir(V*c(\i) +U*(X2)))  -  \r&  A2  =  0,  (35) 


where  T* ,  V ^  and  TZ*  are  adjoint  operators  determined  by  (see  Appendix  A) 


V( \i)=  /  A1(t)M1(C(t),C(s))-Ai(t)M2(C(S), 


dt, 


via j  =  j  (\at)  Maas),at))  -  \i(t)  Maas),at)) 

^(A2)  =  f  ( ^2W 

^  mJo  \ 


dt. 


CM  -  CW 


CW  +  C(«) 


dt. 


Using  the  adjoint  equations-based  method,  we  can  write  the  total  variation  of  the  Lagrangian  only 
with  respect  to  (  and  £: 

SUL  =  Re  |  — (Gq0)  4'C,  r)  -  ( G S°}  C»  +  (Tc(G<0)$C)  +  ^cWC),  AT)  +  <^f(G2<5C),  A2) 

+  (<5cTc(Gf  }C)  +  5(vaG2C),  AT)  +  (  f  ^2  +  ^(G2C),  A2)} 

subject  to  the  constraints  in  (34)  and  adjoint  equations  (35).  Here  we  have 

»i 

'o 


5(r((G^o  =  [  (Gf\t)dt)dMaas),at))  -  g^wcw^cooww)  )  ^ 


<yc^c(C?2C)  =  f1  ( G2(t)C(t )  6M3(as)X(t))  -  G2(t)dt)  <5M4(C(s),  CW))  dt  +  Vc(  ±  tf2  C  «5r), 
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1 


S<;Tk(G 2C)  =-  ,  , 

m  J  0  V 


*r« 


G2(i)C(t)<5 

+  ^c(^2C^r), 


'^(COO.CW)' 

cw  -  as) 


+  G2(i)C(i)  <5 


rn{l\t(s),gt)) 

C(t)  +  C(s) 


dt 


where 


SMj(c t) = *. + Sri + 5z + ?M§i2i  fcl ,  j = 1, 2. 


9r 


Sri 


dz\ 


The  corresponding  derivatives  of  Mj  are  presented  in  Appendix  C. 
We  represented  £(s)  in  the  form  similar  to  (31) 


COO  =  W  X)fe=l  afc  e  WlS  T2fc-2(2s),  s  G  [0,  \  ], 
C(s)  =  C(1  -  «),  se[|,i], 


(36) 


which  in  contrast  to  (31)  allows  for  non-smoothness  at  s  =  | ,  and  calculated  the  gradient  as  in  Section 
3.  Since  the  state  and  adjoint  variables  satisfy  the  following  symmetry  conditions 


H1(l-s)  =  -H1(s),  H2(1-s)  =  H2(s),  H3(1-s)  =  -H3{s), 

Ai(1  -  s)  =  Ai(s),  A2(l  -  s)  =  A2(s), 


we  represented  them  on  [0,  in  the  form  of  series  with  Chebyshev’s  polynomials  of  the  first  kind 
(similar  to  (28)  with  (32))  and  found  unknown  coefficients  by  minimizing  the  total  squared  error  (see 
the  problem  (29)). 

In  our  preliminary  numerical  experiments,  we  obtained  the  optimal  drag  of  0.9877,  normalized  to 
the  drag  of  the  unit- volume  sphere,  for  m  =  9  and  n  =  4  and  expect  that  this  result  can  be  improved. 
The  shape  that  corresponds  to  this  value  is  determined  by 


7  =  0.6192260666309106, 

cti  =  0.9994918803923275,  a2  =  0.18824077143133067, 
a3  =  -0.01408994670305956,  a4  =  0.0066332158943282656, 


and  is  shown  in  Figure  5. 


Figure  5:  Preliminary  results:  optimal  shape  for  the  solid  unit-volume  body  of  revolution  translating 
in  the  fluid  with  constant  velocity  in  the  horizontal  direction  (vertical  axis  is  the  axis  of  revolution, 
and  m  =  9  and  n  =  4). 
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6  Conclusions 


We  have  suggested  the  semi-analytical  approach  to  three-dimensional  (3-D)  shape  optimization  prob¬ 
lems.  The  approach  couples  the  framework  of  the  generalized  analytic  functions  with  adjoint  equations- 
based  method  and  is  summarized  below. 

(i)  Identify  the  class  of  generalized  analytic  functions  related  to  the  governing  equations  (e.g.,  Stokes 
equations,  Oseen  equations,  Maxwell’s  equations,  etc.), 

(ii)  Represent  solution  to  the  governing  equations  in  terms  of  the  generalized  analytic  functions, 

(iii)  Use  the  Cauchy  integral  formula  for  the  generalized  analytic  functions  to  reduce  boundary-value 
problems  for  governing  equations  to  corresponding  integral  equations  (state  equations), 

(iv)  Formulate  shape  optimization  problem  with  integral  equation  constraints  and  apply  adjoint 
equations-based  method  for  obtaining  the  gradient  of  the  Lagrangian, 

(v)  Solve  the  state  and  adjoint  equations  by  minimizing  the  total  squared  error. 

We  have  illustrated  this  approach  in  solving  drag  minimization  problem  for  a  solid  unit-volume 
body  translating  in  a  viscous  incompressible  fluid  with  constant  velocity  under  the  assumption  of  zero 
(low)  Reynolds  number.  We  have  also  solved  this  problem  under  the  constraint  that  the  body  has  the 
axis  of  revolution  and  translates  in  the  direction  perpendicular  to  the  axis.  The  novelty  and  advantage 
of  the  approach  is  in  its  efficiency  and  accuracy,  which  can  be  attributed  to 

■  Solution  on  the  boundary  versus  solution  in  the  domain  (PDE  constrained  optimization), 

■  Analytical  representations  for  the  state  and  adjoint  functions  and  for  body’s  shape, 

■  Analytical  form  of  the  gradient. 

Among  open  issues  are  those  inherited  from  the  gradient-based  method,  namely 

■  Finding  global  optimum  vs.  local  optimum,  and 

■  Finding  initial  shape  efficiently. 

Addressing  these  issues  as  well  as  applying  the  developed  approach  to  other  physical  models,  e.g.,  to 
Maxwell’s  equations  governing  the  behavior  of  electromagnetic  waves,  are  the  subject  for  the  future 
research. 


A  Derivation  of  the  Adjoint  Operator 

Let  <%(/)  =  Jo  (/(*)  Ki(C(s),C(t))  -  fit )  K2(C{s),C(t))Sj  dt  then 


Re{(Sc(/),A>}  =  Re{j\(s)j^  (f(t)  K^s),  C(t))  -  f(t)  K2(<Z(s),  C(*)))  dtdsj 

=  Re{  C  f{s)  f1  \(t)  KiidtfXis))  dtds  -  C  J{s)  P  A(t) 

U  o  Jo  Jo  Jo 

=  Re{Io  So  ^i(C(05  C(s))  —  A(i)  ^(CC^),  C(s))) 

=  Re{(/,Sc*(A))}, 


C(s))  dt  d-s 


and  consequently,  <S£( A)  =  fQL  (\(t)  Ki(((t),  ChO)  -  A(t)  K2(((t),  C(s)))  dt. 
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B  Derivatives  of  Ki(C,  r)  and  /12(C)  r) 


This  section  presents  the  derivatives  of  the  kernels  Ki(^r)  and  i^2(C?r)  entering  the  gradient  expres¬ 
sion. 

nf  K.  r)  +  r?72,r/7l)2  c>«-  T>  n'0,«-  7 ' 

2r  |C  +  t\  J 

me, r)  +  r27?^7l)2  C.K, r) r>)  , 

2ri|C  +  r|  / 

rf°)(C,r)  +  ^-;2,77/l)2C2(C,T)nf«,r))  , 

2r|C-r|  J 

ii!°)(C,  r)  +  r2  m  t 7l)2  e2(C,  T)  Sif  «,7  , 


9Xi(c,7 

_  1  C 

Zl  —  z 

<7(771 

dr 

7ri  C 

l_2(r  -C)2 

2  r 

9Xi(C,7 

_  1  c 

Z  —  Zl 

7(771 

9ri 

7ri  C 

l_2(r  -C)2 

2n 

9X2(C,r ) 

_  1  C 

Z  —  Zl 

7(771 

dr 

7ri  C 

U(c  +  72 

2r 

9X2  (C,  7 

_  1  Cl 

z  -  Zl 

,  7(771 

9ri 


^  VL2(C  +  t) 


2ri 


9/7(77  _  1C  r~ri  .nW(C?r)_ 


9z 


vr;  \2(r  -  C)2 


IC  +  T 


2n  |C  -  t | 

2Ci(77^u;(C,7), 


9//2(77  _  1  /  r  +  n  o(0),  ,  *-*1  r  ((.  _no(0)^ 

g^777  _  9Xi(c,t)  9/7(77  _  dx2(77 

9zi  9z  ’  9zi  dz 

C  Derivatives  of  the  Kernels  for  Transversal  Motion 

This  section  presents  the  derivatives  of  the  kernels  Mi (7  7>  JV(7  7>  M3  (7 7  and  Mi(7  r)  entering 
the  gradient  expression  for  transversal  translation  of  the  body  of  revolution. 


9Mi(77  _  1  C  ~ 

dr 

9Mi(C,7 


Tri  V(r  -  C)2 

1  /  Z  -  Z] 


9ri 


7TV  V(T  -  O2 
dz 


-/c„(C,r))  nfK.rl  +  Z 
+  r)^  )(C>T)  + 


1  r2-r2  +  (z-zi)2^  ^  _,o(0) 


Tri  2r|C  +  r|2 

1  r2  —  r2  +  (z  —  Z1)2 
vri  2ri|C+r|2 


Cn(77^(C,7, 

c,11(77n-)(77. 


Trv(r-C)2  VTX  |C  +  r|2 

9Mi(77  9a/i  (77 


9zi 


9z 


9m2(77 

dr 

9a/2(T7 

9ri 


1  /  ^  -  Ac21(C,r)V°,(C.r)  +  la  +  mjH  e21(C,  r)«f(C.T). 


7ri  \  (r  +  C)2  2r 


Tri  2r  |C  —  t\ 


1/2-21  +  7,(C.T)(  ilW(C.r)  +  7  -7  +  («  ;7>2c21(  <.r)tl!°>(C.r), 


Tri^T  +  C)2  2?’i 


vri  2?’i  |C  —  r| 


9A/2(C,7  _  1  r  +  n  o(0)  _  1  z-zi  o(0)  _ 

- “  “^(FTc?  “  (C’T)  “  SC2l(c’T) jfT7j5n+  (C-T)' 


9z 
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<9M2(C,t)  8M2{(,t) 


dz\ 


dz 


8M3((,t)  _  1  frj-r2 +  (z- Zl)2  -Sn{y{(,T)  r  -  C  -  2r  /o(0) ,  n(i),> 

dr  _  7ri  y  r-C  2r|C  +  r|2  2r(r  -  C)2  V"+  T)  U+^T>)) 

8M3(C,t)  _1  /r2  -r?  +  (z  -  Zl)2  nW (( ,  t)  -  3QW (( ,  t)  t-  (  -  2n  (n(o)fr  ,  0(i),.  +  \ 

dn  7ri  y  r-C  2n|C  +  r|2  +2  n(r  -  C)2  (C’ T)  U+U’t>)) 

8M3(C,t)  _i  Zl-Z  nL0)(C,'T)-3fiL1)(C1r)  + 1 

7T 


7ri  r  —  C 


IC  +  tI2 

8M3((,t)  8M3((,t) 


(+-  o2 


8z\  dz 

8M4((,t)  _1  f  r2  —  r2  +  (z  —  zi)2  n{+\(,T)  +  3g^(CVr)  r  +  C  +  2r  /  (0) 


dr 


m 


1C  -  T\ 


2  r(r  +  C) 


2r(r  +  C)2 


8M4(C,t)  _1^  f  r 2  —  r2  +  (z  —  zi)2  n+\(,r)  +  +  r  +  (  -  2n  fo(0) 


(9?’i  7ri 


1C  -  r| 


2n(r  +  C) 


2ri(r  +  C)2 


(^0)(C,r)  +  ^1}(C,r)) 

(^(C.rJ  +  fi^CC.r) 


8M4((,  r)  _  1  zi-z  nffcr)  +  3fiV'M  _  1  nw(C,r)  +  rJ(C,r) 

7 T 


<9z 


8  | 

^V1}(C,r) 

dr  ' 

V  T"C 

9  | 

(  m+Hm 

<9r3  ' 

{  r-c 

5  ( 

dz  \ 

8  | 

C,r) 

<9r  ' 

V  T"C 

8  | 

^}(C,r) 

<9ri  ' 

^  T-c 

^  / 

dz  \ 

^IC-rl 


1 


T  -  C 

1 

T  -  C 


r  +  C 

c>M4(C,  t)  =  8M4((,t) 

8z\  dz 

or2-r2  +  (z-z  i)2o(i) 


(r  +  C)2 


2r|C  +  T|» 

T  —  C  —  2r-i  (1) 


r  _ri  +  (z_Zl)  T-C-2n  (d  \ 

d  2n|C  •  r|2  -  (C,Tj+  2n(r-C)  +  K’  7  ’ 


r-C  V  IC  +  t 
>(i) 


3^^(C,r)  + 


2ri(r  —  C) 

T~  C  / 


5  /o+J(C>r)\  9  /ri^!}(c,r) 


<9zi  V  r  —  c 

2 

3r'  “  r 


r  +  C 
1 


2r|C  —  r  |2 


3*  V  r-C  )  ’ 


2r(r  +  C) 


+  r)  +  C±4=4l^1,(C,  .) ) , 


r  +  C  V  2ri|C-r|2 


2ri(r +  C) 


—  A*  oMfC  T)  -  in-)(C’r)< 

r  +  C  l  |C-r|2  +  (C’  j  r  +  C 


5  /^a)(C,T)\  9  fny’((,T) 


(1), 


<9+  l  r-C 


3z  l  r-C 
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